\documentclass[twoside,a4paper]{article}
\usepackage{geometry}
\geometry{margin=1.5cm, vmargin={0pt,1cm}}
\setlength{\topmargin}{-1cm}
\setlength{\paperheight}{29.7cm}
\setlength{\textheight}{25.3cm}

% useful packages.
\usepackage[UTF8]{ctex}	% 中文支持
\usepackage{listings} 	% 插入代码
\usepackage{xcolor} 	% 语法高亮
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{enumerate}
\usepackage{graphicx}
\usepackage{multicol}
\usepackage{fancyhdr}
\usepackage{layout}

% some common command
\newcommand{\dif}{\mathrm{d}}
\newcommand{\avg}[1]{\left\langle #1 \right\rangle}
\newcommand{\difFrac}[2]{\frac{\dif #1}{\dif #2}}
\newcommand{\pdfFrac}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\OFL}{\mathrm{OFL}}
\newcommand{\UFL}{\mathrm{UFL}}
\newcommand{\fl}{\mathrm{fl}}
\newcommand{\op}{\odot}
\newcommand{\Eabs}{E_{\mathrm{abs}}}
\newcommand{\Erel}{E_{\mathrm{rel}}}

\newcommand{\R}{\mathbb{R}}
\renewcommand{\proofname}{Proof.}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}{Lemma}[section]
\newtheorem{exercise}{Exercise}[section]
\newtheorem{definition}{Definition}[section]

\lstset{
	columns=fixed,
	numbers=left,                                        % 在左侧显示行号
	numberstyle=\tiny\color{gray},                       % 设定行号格式
	basicstyle=\small,
	tabsize=4,
	frame=single,                                        % 背景边框
	rulecolor=\color[rgb]{0.8,0.8,0.8},                  % 边框颜色
	extendedchars=false,                                 % 解决换页问题
	backgroundcolor=\color[RGB]{245,245,244},            % 设定背景颜色
	keywordstyle=\color[RGB]{40,40,255},                 % 设定关键字颜色
	identifierstyle=\color{black},                       % 普通标识符颜色
	commentstyle=\it\color[RGB]{0,96,96},                % 设置代码注释的格式
	stringstyle=\rmfamily\slshape\color[RGB]{128,0,0},   % 设置字符串格式
	showstringspaces=false,                              % 不显示字符串中的空格
	language=matlab,                                     % 设置语言
	escapeinside=``,                                     % 逃逸字符
	xleftmargin=1em,                                     % 边距设置
	xrightmargin=1em,
	aboveskip=1em
}

\begin{document}
	
	\pagestyle{fancy}
	\fancyhead{}
	\lhead{行一凡 (3190105815)}
	\chead{Project 2}
	\rhead{\today}
	
	\section{项目结构}
	所有程序代码都位于主目录下；三个文件夹 mat、outputs、pictures 依次存放 matlab 绘图脚本、程序输出结果、生成图片。进入主目录下
	\begin{itemize}
		\item 输入 make run 命令编译运行一维程序 oneDimension.cpp 和二维程序 twoDimension.cpp ，输出三个函数样例在多重网格下，每一次 V 循环的收敛速率、总的迭代次数、最终的范数误差、误差的收敛速率，以及 CPU 时间；
		\item 输入 make oneDraw 命令编译运行一维绘图程序 oneDraw.cpp ，输出绘图代码 oneDraw.m ；
		\item 输入 make twoDraw 命令编译运行二维绘图程序 twoDraw.cpp ，输出绘图代码 twoDraw.m ；
		\item 输入 make reduce 命令编译运行测试相对误差的绘图程序 reduce.cpp ，输出绘图代码 reduce.m ；
		\item 输入 make clean 清除可执行程序。
	\end{itemize}
	
	项目使用 JSON for Modern C++ 库读取 json 文件，所需头文件为 json.hpp ，在目录中给出；线性方程组用 LAPACK 求解，通过 lapackSolver.h 头文件调用。
	
	\section{数值实现}
	\subsection{线性系统}
	使用 dim 作为模板参数构造离散线性系统 $ A_{2D}U_{2D}=F_{2D} $，记维数为 $ d $ 。从原点出发，按第一维扫描点，同时检查点在网格中的位置。
	
	如果点不处于边界，则该点对应的矩阵系数设为 $ \frac{2d}{h^2} $ ， $ F_{2D} $ 对应的值设为 $ f(P) $ ，然后对每一个方向，将相邻两点对应的矩阵系数设为 $ -\frac{1}{h^2} $ ，这就完成了标准格点；
	
	如果该点处于边界，但是不处于顶点位置，则使用 ghost cell 近似。按照对应边界的边值条件，不妨设
	\begin{equation}\label{key}
		\alpha u_P + \beta \pdfFrac{u}{n}\bigg|_P = \sigma
	\end{equation}
	将该点对应的矩阵系数设为 $ \alpha+\beta\frac{d}{h} $ ， $ F_{2D} $ 对应值为 $ \sigma +  \beta\frac{h}{2} f(P)  $
	
	如果该点是顶点位置，则使用三点近似
	\begin{equation}\label{key}
		\dfrac{1}{h}\left(\dfrac{3}{2}U_P -2U_Q+\dfrac{1}{2}U_R \right) = \pdfFrac{u}{n}\bigg|_P
	\end{equation}
	其中 $ Q,R $ 是 $ P $ 处法方向反向的相邻两点。此时 $ F_{2D} $ 对应的值设为 $ \sigma $ 。
	 
	\subsection{限制和插值}
	对于 injection 算子，只需要取所有的坐标为偶数的点直接复制即可；对于 full-weighting 算子，按照类似于线性系统的构造方式，从原点出发，按第一维扫描点，同时对该维上的点使用一维 full-weighting 算子。完成扫描后，重新从原点出发，按第二维扫描并对该维上的点使用一维 full-weighting 算子，以此类推即可完成。
	
	对于 linear 算子，类似对每一维进行完全扫描插值即可；最后是 quatratic 算子，利用二次多项式插值，容易得到对 $ x_0,x_1,x_2 $ 点插值 $ f(x_0),f(x_1),f(x_2) $ 后，有
	\begin{equation}\label{key}
		\begin{aligned}
			f\left(\dfrac{x_0+x_1}{2} \right) = \dfrac{1}{8}(3 f(x_0) + 6 f(x_1) - f(x_2))\\
			f\left(\dfrac{x_1+x_2}{2} \right) = \dfrac{1}{8}(- f(x_0) + 6 f(x_1) + 3 f(x_2))
		\end{aligned}
	\end{equation}
	扫描过程类似，插值过程只需要每次取相邻的三个格点代入运算即可。
	
	\section{测试内容}
	\section{测试函数}
	oneDimension.cpp 对以下三个函数
	\begin{equation}\label{key}
		\begin{aligned}
			u(x) &= e^x\\
			u(x) &= \sin x\\
			u(x) &= \cos x
		\end{aligned}
	\end{equation}
	在 $ (0,1) $ 区间上求解边值问题；twoDimension.cpp 对以下三个函数
	\begin{equation}\label{key}
		\begin{aligned}
			u(x) &= e^{y+\sin x}\\
			u(x) &= \sin(x+y) - xy\\
			u(x) &= \sin x + \sin y + xy
		\end{aligned}
	\end{equation}
	在 $ (0,1)^2 $ 区域上求解边值问题。依次使用 Dirichlet、Neumann、Mixed 边值条件，边界数据、使用的限制和插值算子类型、前后松弛的次数、最大迭代次数，以及退出所需的相对误差存放在 inputs.json 文件中。

	\section{数值结果}
	多种数值选项组合的结果存放在 oneRes 和 twoRes 文件中，下面分别展示一维和二维的第一个测试函数的数值结果。使用 V 循环，设置 $ \nu_1=\nu_2=4 $ ，最大迭代次数为 $ 100 $ ，分别测试  Dirichlet、Neumann、Mixed 边值条件，算子组合为 linear + full-weighting 。
	
	\subsection{一维网格}
	不同边值条件下，求解得到的函数图像以及每次 V 循环后残量的收敛图像如下
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.7\linewidth]{../pictures/oneDimension/oneRSInd1}
		\caption{Dirichlet}
		\label{fig:onersind1}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.8\linewidth]{../pictures/oneDimension/oneRSInd2}
		\caption{Neumann}
		\label{fig:onersind2}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.8\linewidth]{../pictures/oneDimension/oneRSInd3}
		\caption{Mixed}
		\label{fig:onersind3}
	\end{figure}

	其中残量图像的纵坐标残量范数取对数，横坐标表示迭代次数。网格结果收敛，并且三种边值条件下迭代次数不超过 20 次。 Dirichlet 边值条件残量收敛最快，约有 4.5 阶收敛速度；其余两种残量收敛较慢，约有 1.5 阶收敛速度。
	
	范数误差收敛图像如下，其中横纵坐标均取对数。三种边界条件下，范数误差都具有 2 阶收敛性。
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.7\linewidth]{../pictures/oneDimension/oneErrorInd1}
		\label{fig:oneerrorind1}
	\end{figure}

	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.7\linewidth]{../pictures/oneDimension/oneTimeInd1}
		\caption{CPU time}
		\label{fig:onetimeind1}
	\end{figure}
	
	逐渐调整相对误差 $ \epsilon $ 减小到 $ 10^{-16} $ ，其中横坐标取对数。可以看出三种边界条件下都不能达到 $ 10^{-16} $ 精度，大约在 $ 10^{-15} $ 到 $ 10^{-13} $ 的相对误差范围。这有可能是限制和插值过程的误差造成的。
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.8\linewidth]{../pictures/reduce/oneEpInd1}
		\caption{Dirichlet}
		\label{fig:oneepind1}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.8\linewidth]{../pictures/reduce/oneEpInd2}
		\caption{Neumann}
		\label{fig:oneepind2}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.7\linewidth]{../pictures/reduce/oneEpInd3}
		\caption{Mixed}
		\label{fig:oneepind3}
	\end{figure}

	\subsection{二维网格}
	不同边值条件下，求解得到的函数图像以及每次 V 循环后残量的收敛图像如下
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.8\linewidth]{../pictures/twoDimension/oneSolutionInd1}
		\caption{Dirichlet}
		\label{fig:onesolutionind1}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.8\linewidth]{../pictures/twoDimension/oneSolutionInd2}
		\caption{Neumann}
		\label{fig:onesolutionind2}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.8\linewidth]{../pictures/twoDimension/oneSolutionInd3}
		\caption{Mixed}
		\label{fig:onesolutionind3}
	\end{figure}
	
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.8\linewidth]{../pictures/twoDimension/oneResidualInd1}
		\caption{Dirichlet}
		\label{fig:oneresidualind1}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.8\linewidth]{../pictures/twoDimension/oneResidualInd2}
		\caption{Neumann}
		\label{fig:oneresidualind2}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.7\linewidth]{../pictures/twoDimension/oneResidualInd3}
		\caption{Mixed}
		\label{fig:oneresidualind3}
	\end{figure}

	其中残量图像的纵坐标残量范数取对数，横坐标表示迭代次数。网格结果收敛，并且三种边值条件下迭代次数不超过 40 次。 Dirichlet 边值条件残量收敛最快，约有 3 阶收敛速度；Neumann 边值条件残量约有 0.1 阶收敛速度； Mixed 边值条件残量约有 1 阶收敛速度。
	
	范数误差收敛图像如下，其中横纵坐标均取对数。三种边界条件下，范数误差都具有 2 阶收敛性。
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.7\linewidth]{../pictures/twoDimension/oneErrorInd1}
		\label{fig:oneerrorind1}
	\end{figure}
	

	逐渐调整相对误差 $ \epsilon $ 减小到 $ 10^{-16} $ ，其中横坐标取对数。可以看出它们都不能达到 $ 10^{-16} $ 精度，大约在 $ 10^{-14} $ 到 $ 10^{-11} $ 的相对误差范围，这有可能是限制和插值过程的误差造成的。特别地对于 Neumann 边值条件，即使迭代 100 次也不能达到 $ 10^{-8} $ 的相对误差，其相对误差在 $ 10^{-3} $ 附近。
	
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.8\linewidth]{../pictures/reduce/twoEpInd1}
		\caption{Dirichlet}
		\label{fig:twoepind1}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.8\linewidth]{../pictures/reduce/twoEpInd3}
		\caption{Mixed}
		\label{fig:twoepind3}
	\end{figure}

	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.7\linewidth]{../pictures/twoDimension/oneTimeInd2}
		\caption{CPU time}
		\label{fig:onetimeind2}
	\end{figure}

	\section{收敛性分析}
	\subsection{线性系统}
	在二维区域 $ \Omega = (0,1)^2 $ 上的 Dirichlet 边值问题
	\begin{equation}\label{key}
		-\Delta u(x) = f(x),\quad u(0) = 0
	\end{equation}
	在网格 $ \Omega^h $ 上离散为 $ m\times m $ 维线性系统
	\begin{equation}\label{key}
		A_{2D}U_{2D} = F_{2D}
	\end{equation}
	其中 $ U_{2D} = (U_{11},U_{21},\cdots,U_{m1},\cdots,U_{1m},U_{2m},\cdots,U_{mm})^T $， $ A_{2D} $ 如下：
	\begin{equation}\label{key}
		A_{2D} = \dfrac{1}{h^2}\left(
		\begin{matrix}
			T & -I\\
			-I & T & -I\\
			& -I & T & -I\\
			& & \ddots & \ddots & \ddots\\
			& & & -I & T & -I\\
			& & & & -I & T
		\end{matrix}
		\right),\quad T = \left(
		\begin{matrix}
			4 & -1\\
			-1 & 4 & -1\\
			& -1 & 4
		\end{matrix}
		\right)
	\end{equation}
	
	\begin{theorem}
		矩阵 $A_{2D}$ 的特征对为
		\begin{equation}\label{key}
			\lambda_{ij} = \lambda_i+\lambda_j,\quad W_{ij} = \mathrm{vec}\left(w_iw_j^T\right),\quad i,j=1,2,\cdots,m
		\end{equation}
		其中 $(\lambda_i,w_i)$ 是 $A$ 的特征对，$ \mathrm{vec} $ 表示 Kronecker 积。
	\end{theorem}
	
	\begin{lemma}
		对于上述线性系统，加权 Jacobi 方法有迭代矩阵
		\begin{equation}\label{key}
			T_{\omega} = (1-\omega)I+\omega D^{-1}(L+U) = I-\dfrac{\omega h^2}{4}A_{2D}
		\end{equation}
		并且它的特征向量与 $ A_{2D} $ 相同，对应的特征值为
		\begin{equation}\label{key}
			\lambda_{ij}(T_{\omega}) = 1-\omega\left(\sin^2\dfrac{i\pi}{2n}+\sin^2\dfrac{j\pi}{2n}\right)
		\end{equation}
	\end{lemma}
	\begin{proof}
		根据分解 $ A_{2D}=D-L-U $ 有
		\begin{equation}\label{key}
			T_\omega = I-\omega D^{-1}(D-L-U) = I-\omega D^{-1}A_{2D} = I-\dfrac{\omega h^2}{4}A_{2D}
		\end{equation}
		其中 $ D^{-1} $ 是对角阵，且对角元为 $\frac{h^2}{4}$ ；代入特征对有
		\begin{equation}\label{key}
			T_\omega W_{ij} = \left(I-\dfrac{\omega h^2}{4}A_{2D}\right)W_{ij} = \left(1-\dfrac{\omega h^2}{4}\lambda_{ij}(A_{2D})\right)W_{ij} = \left(1-\omega\left(\sin^2\dfrac{i\pi}{2n}+\sin^2\dfrac{j\pi}{2n}\right)\right)W_{ij}
		\end{equation}
		对于 $ i,j=1,2,\cdots,n-1 $ 成立，即证。
	\end{proof}
	
	\subsection{限制和插值}
	\begin{definition}
		限制算子
		\begin{equation}\label{key}
			I_h^{2h}:V^h\mapsto V^{2h}
		\end{equation}
		将细网格 $ \Omega^h $ 映射到粗网格 $ \Omega^{2h} $ ，其中 $ V^h $ 表示二维网格点。
	\end{definition}
	
	\begin{definition}
		injection 算子是限制算子
		\begin{equation}\label{key}
			V_{ij}^{2h} = V_{2i,2j}^h 
		\end{equation}
		其中 $ i,j=1,\cdots,\frac{n}{2}-1 $ 。
	\end{definition}
	
	\begin{definition}
		full-weighting 算子是限制算子
		\begin{equation}\label{key}
			\begin{aligned}
				V_{ij}^{2h} = \dfrac{1}{16}&(V_{2i-1,2j-1}^h + 2V_{2i,2j-1}^h + V_{2i+1,2j-1}^h\\
				&+ 2V_{2i-1,2j}^h + 4V_{2i,2j}^h + 2V_{2i_1,2j}^h\\
				& + V_{2i-1,2j+1}^h + 2V_{2i,2j+1}^h + V_{2i+1,2j+1}^h)
			\end{aligned}
		\end{equation}
		其中 $ i,j=1,\cdots,\frac{n}{2}-1 $ 。
	\end{definition}
	
	\begin{definition}
		插值算子
		\begin{equation}\label{key}
			I_{2h}^{h}:V^{2h}\mapsto V^{h}
		\end{equation}
		将粗网格 $ \Omega^{2h} $ 映射到细网格 $ \Omega^{h} $ ，其中 $ V^h $ 表示二维网格点。
	\end{definition}
	
	\begin{definition}
		linear interpolation 算子是插值算子
		\begin{equation}\label{key}
			\begin{aligned}
				V_{2i,2j}^{h} &= V_{ij}^{2h}\\
				V_{2i+1,2j}^h &= \frac{1}{2}(V_{i,2j}^{2h}+V_{i+1,2j}^{2h})\\
				V_{2i,2j+1}^h &= \frac{1}{2}(V_{2i,j}^{2h}+V_{2i,j+1}^{2h})\\
				V_{2i+1,2j+1}^h &= \frac{1}{4}(V_{i,j}^{2h}+V_{i+1,j}^{2h}+V_{i,j+1}^{2h}+V_{i+1,j+1}^{2h})
			\end{aligned}
		\end{equation}
		其中 $ i,j=1,\cdots,\frac{n}{2}-1 $ 。
	\end{definition}
	
	\subsection{收敛性}
	\begin{lemma}
		互补模式满足
		\begin{equation}\label{key}
			\begin{aligned}
				(W_{i^\prime j}^h)_{(p,q)} &= (-1)^{p+1}(W_{ij}^h)_{(p,q)}\\
				(W_{i j^\prime}^h)_{(p,q)} &= (-1)^{q+1}(W_{ij}^h)_{(p,q)}\\
				(W_{i^\prime j^\prime}^h)_{(p,q)} &= (-1)^{p+q}(W_{ij}^h)_{(p,q)}\\
			\end{aligned}
		\end{equation}
	\end{lemma}
	\begin{proof}
		利用 $ (W_{ij}^h)_{(p,q)} = w_{i,p}^hw_{j,q}^h $ 以及 Lemma 9.37 代入即证。
	\end{proof}
	
	\begin{lemma}
		full-weighting 算子作用于 $ \Omega^h $ 上的互补模式有
		\begin{equation}\label{key}
			\begin{aligned}
				I_{h}^{2h}W_{ij}^h &= c_ic_jW_{ij}^{2h} = \cos^2\dfrac{i\pi}{2n}w_i^{2h}\cdot\cos^2\dfrac{j\pi}{2n}w_j^{2h}\\
				I_{h}^{2h}W_{i^\prime j}^h &= -s_ic_jW_{ij}^{2h} = -\sin^2\dfrac{i\pi}{2n}w_i^{2h}\cdot\cos^2\dfrac{j\pi}{2n}w_j^{2h}\\
				I_{h}^{2h}W_{ij^\prime}^h &= -c_is_jW_{ij}^{2h} = -\cos^2\dfrac{i\pi}{2n}w_i^{2h}\cdot\sin^2\dfrac{j\pi}{2n}w_j^{2h}\\
				I_{h}^{2h}W_{i^\prime j^\prime}^h &= s_is_jW_{ij}^{2h} = \sin^2\dfrac{i\pi}{2n}w_i^{2h}\cdot\sin^2\dfrac{j\pi}{2n}w_j^{2h}
			\end{aligned}
		\end{equation}
		其中 $ i,j\in[1,\frac{n}{2}),\ i^\prime=n-i,\ j^\prime=n-j $ ，以及
		\begin{equation}\label{key}
			\forall k\in\left[1,\frac{n}{2}\right),\quad I_h^{2h}W_{k,\frac{n}{2}}^h = I_h^{2h}W_{\frac{n}{2},k}^h = 0
		\end{equation}
	\end{lemma}
	\begin{proof}
		利用定义，我们有
		\begin{equation}\label{key}
			\begin{aligned}
				(I_h^{2h}W_{ij}^h)_{(p,q)} &= \dfrac{1}{16}\bigg(\sin\left(\dfrac{2p-1}{n}i\pi\right)\sin\left(\dfrac{2q-1}{n}j\pi\right) + 2\sin\left(\dfrac{2p}{n}i\pi\right)\sin\left(\dfrac{2q-1}{n}j\pi\right) + \sin\left(\dfrac{2p+1}{n}i\pi\right)\sin\left(\dfrac{2q-1}{n}j\pi\right)\\
				&+ 2\sin\left(\dfrac{2p-1}{n}i\pi\right)\sin\left(\dfrac{2q}{n}j\pi\right) + 4\sin\left(\dfrac{2p}{n}i\pi\right)\sin\left(\dfrac{2q}{n}j\pi\right) + 2\sin\left(\dfrac{2p+1}{n}i\pi\right)\sin\left(\dfrac{2q}{n}j\pi\right)\\
				&+ \sin\left(\dfrac{2p-1}{n}i\pi\right)\sin\left(\dfrac{2q+1}{n}j\pi\right) + 2\sin\left(\dfrac{2p}{n}i\pi\right)\sin\left(\dfrac{2q+1}{n}j\pi\right) + \sin\left(\dfrac{2p+1}{n}i\pi\right)\sin\left(\dfrac{2q+1}{n}j\pi\right) \bigg)\\
				&= \dfrac{1}{16}\left(2+2\cos\dfrac{i\pi}{n}\right)\sin\left(\dfrac{2p}{n}i\pi\right)\bigg(\sin\left(\dfrac{2q-1}{n}j\pi\right) + 2\sin\left(\dfrac{2q}{n}j\pi\right) + \sin\left(\dfrac{2q+1}{n}j\pi\right)\bigg)\\
				&= \dfrac{1}{16}\left(2+2\cos\dfrac{i\pi}{n}\right)\sin\left(\dfrac{2p}{n}i\pi\right)\left(2+2\cos\dfrac{j\pi}{n}\right)\sin\left(\dfrac{2q}{n}j\pi\right)\\
				&= \cos^2\dfrac{i\pi}{2n}w_i^{2h}\cdot\cos^2\dfrac{j\pi}{2n}w_j^{2h}
			\end{aligned}
		\end{equation}
		这就证明了第一个式子，其余三个等式进行替换 $ i\to n-i,\ j\to n-j $ 后类似地得到。
	\end{proof}
	
	\begin{lemma}
		linear interpolation 算子作用于 $ \Omega^{2h} $ 有
		\begin{equation}\label{key}
			I_{2h}^hW_{ij}^{2h} = c_ic_jW_{ij}^h-c_is_jW_{ij^\prime}^h-s_ic_jW_{i^\prime j}^h+s_is_jW_{i^\prime j^\prime}^h
		\end{equation}
	\end{lemma}
	\begin{proof}
		利用互补模式的性质，将上面右端整理得到
		\begin{equation}\label{key}
			RHS_{(p,q)} = \left(\cos^2\dfrac{i\pi}{2n}+(-1)^p\sin^2\dfrac{i\pi}{2n} \right)\left(\cos^2\dfrac{j\pi}{2n}+(-1)^q\sin^2\dfrac{j\pi}{2n} \right)w_{i,p}^hw_{j,q}^h
		\end{equation}
		另一方面根据定义有
		\begin{equation}\label{key}
			(I_h^{2h}W_{ij}^h)_{(p,q)} = \left\{
			\begin{matrix}
				w_{i,p}^hw_{j,q}^h & p,q \mid 2\\
				\cos\frac{i\pi}{n}w_{i,p}^hw_{j,q}^h & p\nmid 2,\ q\mid 2\\
				\cos\frac{j\pi}{n}w_{i,p}^hw_{j,q}^h & p\mid 2,\ q\nmid 2\\
				\cos\frac{i\pi}{n}w_{i,p}^h\cdot\cos\frac{j\pi}{n}w_{j,q}^h & p,q\nmid 2
			\end{matrix}
			\right.
		\end{equation}
		容易看出，按照上面 $ p,q $ 的不同情况，得到的 $ RHS $ 与 $ (I_h^{2h}W_{ij}^h)_{(p,q)} $ 分别对应，即证。
	\end{proof}
	
	\begin{theorem}
		双网格校正算子在子空间 $ \mathbb{W}_{ij}^h = \mathrm{span} \{W_{ij}^h,W_{i^{\prime}j}^h,W_{ij^{\prime}}^h,W_{i^{\prime}j^{\prime}}^h\} $ 上是不变算子，满足
		\begin{equation}\label{key}
			\begin{aligned}
				TGW_{ij} &= \lambda_{ij}^{\nu_1+\nu_2}(1-c_ic_jK_1)W_{ij}^h+
				\lambda_{ij}^{\nu_1}\lambda_{ij^\prime}^{\nu_2}c_is_jK_1W_{ij^\prime}^h+
				\lambda_{ij}^{\nu_1}\lambda_{i^\prime j}^{\nu_2}s_ic_jK_1W_{i^\prime j}^h-
				\lambda_{ij}^{\nu_1}\lambda_{i^\prime j^\prime}^{\nu_2}s_is_jK_1W_{i^\prime j^\prime}^h\\
				TGW_{ij^\prime} &= \lambda_{ij^\prime}^{\nu_1}\lambda_{i j}^{\nu_2}c_ic_jK_2W_{ij}^h+
				\lambda_{ij^\prime}^{\nu_1+\nu_2}(1-c_is_jK_2)W_{ij^\prime}^h-
				\lambda_{ij^\prime}^{\nu_1}\lambda_{i^\prime j}^{\nu_2}s_ic_jK_2W_{i^\prime j}^h+
				\lambda_{ij^\prime}^{\nu_1}\lambda_{i^\prime j^\prime}^{\nu_2}s_is_jK_2W_{i^\prime j^\prime}^h\\
				TGW_{i^\prime j} &= \lambda_{i^\prime j}^{\nu_1}\lambda_{i j}^{\nu_2}c_ic_jK_3W_{ij}^h-
				\lambda_{i^\prime j}^{\nu_1}\lambda_{ij^\prime}^{\nu_2}c_is_jK_3W_{ij^\prime}^h+
				\lambda_{i^\prime j}^{\nu_1+\nu_2}(1-s_ic_jK_3)W_{i^\prime j}^h+
				\lambda_{i^\prime j}^{\nu_1}\lambda_{i^\prime j^\prime}^{\nu_2}s_is_jK_3W_{i^\prime j^\prime}^h\\
				TGW_{i^\prime j^\prime} &= \lambda_{i^\prime j^\prime}^{\nu_1}\lambda_{i j}^{\nu_2}c_ic_jK_4W_{ij}^h+
				\lambda_{i^\prime j^\prime}^{\nu_1}\lambda_{ij^\prime}^{\nu_2}c_is_jK_4W_{ij^\prime}^h+
				\lambda_{i^\prime j^\prime}^{\nu_1}\lambda_{i^\prime j}^{\nu_2}s_ic_jK_4W_{i^\prime j}^h-
				\lambda_{i^\prime j^\prime}^{\nu_1+\nu_2}(1-s_is_jK_4)W_{i^\prime j^\prime}^h
			\end{aligned}
		\end{equation}
		其中
		\begin{equation}\label{key}
			K_1 = \dfrac{c_ic_j(s_i+s_j)}{s_ic_i+s_jc_j},\ 
			K_2 = \dfrac{s_ic_j(c_i+s_j)}{s_ic_i+s_jc_j},\ 
			K_3 = \dfrac{c_is_j(s_i+c_j)}{s_ic_i+s_jc_j},\ 
			K_4 = \dfrac{s_is_j(c_i+c_j)}{s_ic_i+s_jc_j}
		\end{equation}
		以及 $ \lambda_{ij} $ 为 $ T_\omega $ 的特征值。
	\end{theorem}
	\begin{proof}
		先考虑 $ \nu_1=\nu_2=0 $ 的情况，就有
		\begin{equation}\label{key}
			\begin{aligned}
				A_{2D}^hW_{ij}^h &= \dfrac{4(s_i+s_j)}{h^2}W_{ij}^h\\
				I_h^{2h}A_{2D}^hW_{ij}^h &= \dfrac{4(s_i+s_j)c_ic_j}{h^2}W_{ij}^{2h}\\
				(A_{2D}^{2h})^{-1}A_{2D}^hW_{ij}^h &= \dfrac{4(s_i+s_j)c_ic_j}{h^2}\dfrac{(2h)^2}{4(\sin^2\frac{i\pi}{n}+\sin^2\frac{j\pi}{n})}W_{ij}^{2h}\\
				-I_{2h}^h(A_{2D}^{2h})^{-1}A_{2D}^hW_{ij}^h &= -\dfrac{(s_i+s_j)c_ic_j}{s_ic_i+s_jc_j}(c_ic_jW_{ij}^h-c_is_jW_{ij^\prime}^h-s_ic_jW_{i^\prime j}^h+s_is_jW_{i^\prime j^\prime}^h)\\
				[I-I_{2h}^h(A_{2D}^{2h})^{-1}A_{2D}^h]W_{ij}^h &= (1-c_ic_jK_1)W_{ij}^h+c_is_jK_1W_{ij^\prime}^h+s_ic_jK_1W_{i^\prime j}^h-s_is_jK_1W_{i^\prime j^\prime}^h
			\end{aligned}
		\end{equation}
		类似地有
		\begin{equation}\label{key}
			\begin{aligned}
				A_{2D}^hW_{i^\prime j}^h &= \dfrac{4(s_{i^\prime}+s_j)}{h^2}W_{i^\prime j}^h = \dfrac{4(c_{i}+s_j)}{h^2}W_{i^\prime j}^h\\
				I_h^{2h}A_{2D}^hW_{i^\prime j}^h &= -\dfrac{4(c_i+s_j)s_ic_j}{h^2}W_{ij}^{2h}\\
				(A_{2D}^{2h})^{-1}A_{2D}^hW_{i^\prime j}^h &= -\dfrac{4(c_i+s_j)s_ic_j}{h^2}\dfrac{(2h)^2}{4(\sin^2\frac{i\pi}{n}+\sin^2\frac{j\pi}{n})}W_{ij}^{2h}\\
				-I_{2h}^h(A_{2D}^{2h})^{-1}A_{2D}^hW_{i^\prime j}^h &= \dfrac{(c_i+s_j)s_ic_j}{s_ic_i+s_jc_j}(c_ic_jW_{ij}^h-c_is_jW_{ij^\prime}^h-s_ic_jW_{i^\prime j}^h+s_is_jW_{i^\prime j^\prime}^h)\\
				[I-I_{2h}^h(A_{2D}^{2h})^{-1}A_{2D}^h]W_{i^\prime j}^h &= c_ic_jK_2W_{ij}^h-c_is_jK_2W_{ij^\prime}^h+(1-s_ic_jK_2)W_{i^\prime j}^h+s_is_jK_2W_{i^\prime j^\prime}^h
			\end{aligned}
		\end{equation}
		其余两式同理可得。
		
		最后考虑 $ T_{\omega} $ 作用，它与 $ A_{2D} $ 有相同特征向量，且有特征值 $ \lambda_{ij} = 1-\omega(s_i+s_j) $ 。不妨先分析第一个等式，$ T_\omega^{\nu_1} $ 预先作用于 $ W_{ij} $ 得到系数 $ \lambda_{ij}^{\nu_1} $ ；之后 $ T_\omega^{\nu_2} $ 作用于上面分析后的结果，因此对每一项 $ W_{ij}^h,W_{i^{\prime}j}^h,W_{ij^{\prime}}^h,W_{i^{\prime}j^{\prime}}^h $ 得到对应其特征值的系数 $ \lambda_{ij}^{\nu_2},\lambda_{i^{\prime}j}^{\nu_2},\lambda_{ij^{\prime}}^{\nu_2},\lambda_{i^{\prime}j^{\prime}}^{\nu_2} $ 。
	\end{proof}
	
	\begin{exercise}
		将上面的等式改写为矩阵形式
		\begin{equation}\label{key}
			TG \begin{pmatrix}
				W_{ij}\\
				W_{ij^\prime}\\
				W_{i^\prime j}\\
				W_{i^\prime j^\prime}
			\end{pmatrix} = 
			\begin{pmatrix}
				c_{11} & c_{12} & c_{13} & c_{14} \\
				c_{21} & c_{22} & c_{23} & c_{24} \\
				c_{31} & c_{32} & c_{33} & c_{34} \\
				c_{41} & c_{42} & c_{43} & c_{44} 
			\end{pmatrix}
			\begin{pmatrix}
				W_{ij}\\
				W_{ij^\prime}\\
				W_{i^\prime j}\\
				W_{i^\prime j^\prime}
			\end{pmatrix}
		\end{equation}
		从前面得到的系数可以看出，每个 $ K_i $ 都不大于 $ 1 $ 。若 $ i,j $ 都对应低频，则矩阵中第一个元素中 $ 1-c_ic_jK_1 $ 很小，剩余元素由于包含 $ s_i,s_j $ 或者 $ \lambda_{i^\prime j},\lambda_{ij^\prime},\lambda_{i^\prime j^\prime} $ 因此都很小，这就保证了收敛性。
	\end{exercise}
	
\end{document}